On the nature of interlayer interactions in a system of two graphene fragments 
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With the help of the quantum chemistry methods we have investigated the nature of interlayer 
interactions between graphene fragments in different stacking arrangements ( AA and AB) . We found 
that the AB stacking pattern as the ground state of the system, is characterized by the effective 
inter-band orbital interactions which are barely present in the AA. Their vanishing induces electronic 
decoupling between the graphene layers, so that the bonding interaction AEoi between the flakes 
is drastically reduced from —0.482 eV to —0.087 eV as the stacking pattern is changed from AB to 
AA. The effective way to improve the bonding interaction between layers preserving the same A A 
lattice order is to induce rotation of the layer. As the flake is rotated, the bonding interactions are 
improved mostly due to suppression of the Pauli repulsion which in turn increases the interlayer 
orbital interactions, while the inter-band part of those remain negligible on the whole range of the 
rotation angle. The Pauli repulsion is also found to be the main force moving apart two fragments as 
the stacking pattern is changed from AA to AB. This enhances the equilibrium interlayer distance, 
which for the AA staking is larger than the established value for the AB stacking (3.4 A). 



three layers, the interlayer spacing depends on the num- 
ber of layers [l5| due to different nodal interactions of the 
overlapping n orbitals. For an odd number of adjacent 
layers the equilibrium spacing between the layers was pre- 
dicted to be 3.30 A, while for the even number it is 3.58 
A [15]. It then clearly follows that in graphite the elec- 
tronic coupling between the layers can not be neglected. 
Rather, it provides a substantial influence on the inter- 
layer interactions in addition to the van der Walls and 
steric type of interactions. 

In this context, it is worth pondering what is acually 
happening with the interlayer interactions in twisted bi- 
layer graphene where interestingly, one observes elec- 
tronic decoupling between the layers. The decoupling 
was first observed experimentally [ll-Q and was later in- 
vestigated theoretically |lll-[l3|. Theoretical interpreta- 
tions |lll - [l3j relate the decoupling to the occurrence of 
a misorientation of 2°-5° between the layers. According 
to a proposed model, the layer rotation in the real space 
induces a displacement of the Dirac cones generated in 
each layer in the reciprocal space |ll| - [l3| thereby caus- 
ing the interlayer decoupling. Experimentally, the mis- 
orientation of 2°-5° in AA-stacked bilayer graphene has 
been detected in systems that were created using vari- 
ous fabrication techniques such as the epitaxial growth 
l-Q, chemical vapor deposition [J| and ultrasonication 
5|. Appearance of the rotational misorientation that is 
independent of the fabrication techniques suggests the 
presence of some forces between the layers strong enough 
to cause the layer rotation. However, the available theo- 
retical models [lll - [l3| on the electronic properties of the 
AA stacked graphene deal only with the band properties 
but do not shed any light on the underlying physical rea- 
sons involved in decoupling, such as the interlayer forces. 
In a recent work [18'] we related the origin of the de- 
coupling phenomenon and rotational misorientation with 
layer stacking pattern which is AA in fabricated multi- 
layer graphene [IrtSi] against the AB stacking in natural 
graphite. For the AA staking, the interlayer electronic 
coupling is suppressed by a significant repulsion arising 



Bilayer graphene [1-|7[ , being regarded as an important 
system for applications in semiconductor electronics has 
intrigued the scientific community to look deeper into the 
nature of interlayer interactions in this system. Graphene 
layers stacked together is uniquely different from other 
solid state materials because of the interlayer weak van 
der Waals and steric interactions, instead of the occur- 
rence of C-C bonding between the layers. The half- filled 
p orbital left on each carbon atom after bonding with 
its neighbors within the honeycomb lattice is responsible 
for formation of the n bonds between two neighboring 
atoms within the same layer. This creates a closed shell 
electron system carrying the weakly bound n electrons 
that are distinct due to their high mobility within the 
graphene layer |8l-[l0|. Stacking two systems of closed 
electron shells would cause the interlayer interactions to 
be mostly repulsive which results in the expulsion of va- 
lence electrons from the overlap region such that only a 
weak electronic coupling between the layers can occur. If 
under certain conditions the electronic coupling becomes 
negligible then each layer would display its own electronic 
behavior in the band diagram as was found recently in 
twisted bilayer graphene |ll| - [l3| . 

Even before the discovery of graphene the interlayer 
interactions in natural graphite, which is basically a sys- 
tem of stacked graphene layers, received intense atten- 
tion J14|4l7| where the effect of interlayer decoupling was 
not encountered. Three types of layer arrangements are 
known to exist in graphite [Tjj , but the most common one 
is the Bernal stacking in which the carbon atoms belong- 
ing to different sublattices A and B form the AB stacking 
pattern between the layers. Contrary to the conventional 
wisdom that only the long-range van der Waals interac- 
tion is important in the case of stacking of closed shell 
systems, for the AB stacking in graphite it was shown 
that the orbital overlap between the n orbitals belonging 
to different layers [1^, [l^ is as essential as the all impor- 
tant van der Waals forces. Therefore, it was predicted 
that despite the well known interlayer distance of 3.4 A 
in natural graphite, for a system containing only two or 



between the graphene layers [18|. This repulsion is also 
expected to be responsible for the occurrence of lattice 
misorientation between the layers. It was suggested that 
rotational misorientation, which creates the Moire pat- 
tern, appears as a way to suppress the repulsion, thereby 
lowering the total energy of the system. Even a slight 
layer rotation of ^2°-5° substantially shrinks the areas 
characterized by the AA lattice superposition in which 
repulsion dominates over other forces (the areas with AA 
and AB stacking coexist in the Mo ire p attern) . The other 
important result was a prediction |l8] that the strong re- 
pulsion may induce bumps on the graphene surface in 
the areas where AA stacking is preserved. All these ef- 
fects are important and require careful studies because 
the phenomenon of layer rotation through electronic cou- 
pling between layers can offer ways to manipulate the 
electronic properties of twisted graphene (Moire pattern 
of different rotation angle is characterized by different 
percentage of AA-spotted areas). Therefore, in this work 
we present a detailed quantitative analysis of the repul- 
sive forces and the orbital overlap in stacked graphene 
layers and their alteration with the appearance of rota- 
tion. Our studies are based on the density functional 
methods including a recently proposed empirical correc- 
tion (Grimme correction [l9|) which was developed for 
a proper consideration of the dispersive interactions be- 
tween the closed shell electron systems. 



I. COMPUTATIONAL METHODS 

The computations were performed with the ADF quan- 
tum chemistry code [20| which uses the Kohn-Sham ap- 
proach to density functional theory (DFT). The Kohn- 
Sham approach replaces the many-body system within 
the Hamiltonian equation by a system of the non- 
interacting particles while all the many body terms are 
incorporated into the so-called Kohn-Sham potential. 
This concept is quite useful in the investigation of inter- 
acting closed shell systems because it allows us to present 
each graphene flake as an isolated fragment and two frag- 
ments interacts as the flakes are stacked. In this way, a 
proper investigation of the forces and the orbital over- 
lap can be performed directly in terms of the fragment 
presentation. 

Within the ADF code the forces between frag- 
ments are included in the bonding energy AE'^ which 
comprises of several majors components [21| (A_E°= 
AVei+AEp+AEprep+AEoi+Edis)- The first component 
(AVei) takes care of the interactions of electrostatic na- 
ture related to the modification of the charge distribution 
(originated from the charge transfer between occupied 
and unoccupied orbitals), when two systems are allowed 
to interact. The second one is the energy change induced 
by the Pauli repulsion (AEp), which include several com- 
ponents; exchange repulsion, kinetic repulsion, overlap 
repulsion, all results from obeying the Pauli antisymme- 
try principle. The next term AEprep describes the energy 



required to change the conformation of the fragments 
(structural modification) from the initial geometry con- 
taining separate fragments to the final geometry where 
the fragments are allowed to interact. The bonding inter- 
actions between two fragments are included in the AEoi 
term which originates from the overlap of the fragment's 
orbitals. The last term Edis is the empirical dispersion 
correction introduced by Grimme I.19l| and its magnitude 
is defined by the long-range van der Waals interactions, 
whose contribution in the short range is reduced by the 
damping function. 

Even though the AEoi term is a measure of the orbital 
overlap, the interlayer forces such as the Pauli repulsion 
and orbital polarization contribute to the AEoi as well. 
The effect of interlayer forces can not be discarded from 
AEoi and so the overlap of the selected orbitals can not 
be separated from the others. This makes it hard to get 
a proper understanding of the intricacies of interlayer in- 
teractions between two fragments. The most effective 
way to proceed is to follow the established method of 
linear combinations of the orbitals. This method can be 
applied for the results obtained with the ADF program. 
With the ADF, each fragment is described by its own set 
of orbitals and the program facilitates their mixing upon 
the inclusion of the interaction between the fragments. 
Therefore, the fragment approach allows us to evaluate 
the overlap matrix Si^j between the fragments i and j 
of the Kohn-Sham Hamiltonian (('/'J^ksI'^;)) directly in 
terms of the linear combinations of the fragment orbitals 
via the relation /iks = SCEC~^, where C is the eigen- 
vector defined in terms of the fragment orbitals and E is 
the eigenvalue matrix [22l |. The overlap matrix S purely 
depends on the form of the interacting orbitals and on the 
distance that keeps the two fragments apart neglecting 
the contribution from the attractive and repulsive forces 
arising between the fragments. We used the overlap ma- 
trix Sij to define the spatial overlap integral between the 
fragments i and j, which is Jij — {lp^\H\lpj). 

In this work we consider the spatial overlap integral 
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between the highest occupied fragment orbitals 
(HOFO), i.e., between two tt orbitals, each located on dif- 
ferent fragments while their overlap defines the HOMO of 
the joint system. The overlap integral was also calculated 
between the it and tt* orbitals, i.e., between the highest 
occupied orbital of one fragment (HOFOi) with the low- 
est unoccupied orbital of another fragment (LUFO2) and 
because there are two parts of such interactions, HOFOi- 
LUFO2 and HOFO2-LUFO1, the average value of overlap 
integral was considered and combined into the J^ ^ . 

We used the hybrid BLYP exchange-correlation func- 
tional, applying the empirical dispersion correction 1.05 
recommended by Grimme [l9|, |23|. For the interacting 
molecules of closed electron shells it was found that the 
proposed correction is enough to reproduce the inter- 
molecular distance to what is observed in the experiments 
or achieved with a more accurate level such as the ah ini- 
tio M0ller-Plesset second-order (MP2) method [11]. For 
a proper description of the tails of the electron wave- 



functions that is important for long-range interactions, 
we used the Slater-type orbitals. The quite extended 
TZP basis set (triple-^ polarized basis set) was applied 
in all the calculations which improves the precision of our 
results while suppressing the basis set superposition er- 
We tested the chosen method to reproduce the 
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interlayer distance between the graphene flakes stacked 
in the AB pattern (which is well known to be 3.4 A in 
natural graphite) and indeed the correct interlayer dis- 
tance was obtained (the so-called equilibrium distance of 
the AB pattern rfegfAS))- For this calculation the atomic 
coordinates within the graphene plane (along the x, y di- 
rections) were confined and only the coordinates in the z 
direction, i.e., perpendicular to the graphene plane, were 
used for relaxation. In fact, the full relaxation of the 
system of two stacked fragments is problematic as the 
repulsion between layers leads to sliding of the fragments 
away from each other. 

Our main emphases in this work are the energy de- 
composition analysis of the bonding energy AE^ {AVei, 
AEp, AEprep and AEoi) and the spatial overlap inte- 
grals {J^~^ , Jilf^) which are determined by the stack- 
ing pattern between two graphene flakes (flake rotation) 
and the interlayer distance. In most cases the single point 
calculations have been used for which the contribution of 
AEprep becomes zero. 



II. AA AND AB STACKING 

For our investigations we used two graphene flakes with 
the carbon atoms at the edges terminated by the hy- 
drogen atoms, as shown in Fig. [T] (a) with our goal to 
minimize the contribution of the localized states into the 
simulation results. Because bilayer graphene obtained in 
the experiments [1-^ has shown the AA stacking pattern 
instead of the AB, common in natural graphite, we probe 
the AA stacking for the equilibrium distance between the 
flakes (rfeg(AA))- It was found that this distance is in- 
deed enhanced in the AA stacking up to deq(yiA)=3.67 
A against the deq(AB)=3.4 A for the AB stacking. To 
explain an increase in the interlayer distance we applied 
the decomposition analysis of the bonding energy. 

We performed the single point calculations in which 
the layers were separated by the established equilibrium 
distances (deg) found to be different for A A and AB 
stacking. The obtained interlayer forces and spatial over- 
lap integrals are collected in Table HI For the purpose of 
comparison, in addition to the lattice arrangements AA 
and AB, we also carried out the bonding analysis for the 
AA' stacking pattern for which the interlayer equilibrium 
distance was found to be 3.41 A and those results are also 
enclosed in Table HI The AB-conformation is character- 
ized by a much stronger bonding interaction than that of 
AA, as defined by a more negative bonding energy AE^ ^ 
thereby making the AB-stacking the ground state of the 
system. The structural distinction of the AB configura- 
tion from AA consists in sliding of one graphene flake 
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FIG. 1: (a) Stacked graphene flakes in the 3D space coordi- 
nate system, (b) The energy components {AVei, AEp,AEoi) 
of the bonding energy (A_B'^) as a function of the interlayer 
distance between two graphene fragments. AVei is the elec- 
trostatic interactions, AEp is Pauli repulsion, AEoi is the 
orbital interactions energy. 



relative to the other by 1.42 A along the x axis that in- 
duces the shift of the bond positions between the layers 
against their superposition for the AA stacking (for the 
AA case the atomic coordinates are matched in the x and 
y directions for both layers). To obtain the AA' stacking 
a sliding of 1.23 A is applied along the y axis. 

For the AA stacking, due to the lattice superposition 
the TT clouds between the layers are also superposed that 
leads to their effective overlap defined by the spatial over- 
lap integral jf-^^—0AA3 eV. The flake sliding induced 
for the AB and AA' stacking breaks the bond superpo- 
sition condition recognized for the AA stacking thereby 
inducing the disarrangement of the n clouds within the 
overlapping region. As a result, the tt — tt interlayer 
interactions, which is of particular interest since it is 
supposedly responsible for the occurrence of electronic 
coupling between the graphene layers [3, [23 , is signifi- 
cantly reduced for the AB stacking pattern. This is re- 
flected by a suppression of the spatial overlap integral to 
J, ~ =—0.229 eV (even more drastic reduction is ob- 

served for the AA' case where j"f" =0.023 eV). 

When two fragments are stacked, the majority of the 
orbital interactions described by AEoi arise from the 
overlap of the tt orbitals (tt-tt or tt-tt*). For interaction 
of two closed shell systems, the tt-tt overlap is not the 



TABLE I: The electronic properties and the interlayer forces 
between two graphene fragments stacked in AA (de(j(AA)=3.67 
A), in AB (de9(AS)=3.4 A) and in AA'{deq(AA)=SAl A) con- 
figurations. All values are in eV. 





HOMO 


LUMO 


jH-H 


tH-L 


AEo^ 


AVe, 


AEp 


AS" 


AA 
AB 

AA' 


-4.024 
-4.164 
-4.070 


-2.596 
-2.506 
-2.565 


0.443 
-0.229 
0.023 


10-" 
-0.169 
-0.191 


-0.087 
-0.482 
-0.491 


-0.691 
-1.184 
-1.263 


2.136 
3.388 
3.644 


-1.483 
-1.931 
-1.863 



one that leads to the bonding interactions, and there- 
fore, might be ignored within the orbital interaction term 
AEoi- That explains the contradictory behavior of the 
overlap integral J^f^ and the AEoi term, such that 

when a reduction of the overlap integral J^f^ occurs for 
the AB stacking, the orbital interactions AEoi between 
the fragments is improved. However, the improvement of 
AEoi with modification of the lattice arrangement from 
AA to AB is also not consistent with behavior of the Pauli 
repulsion AEp whose increase supposedly suppresses the 
interlayer orbital interaction AEoi- Therefore, to under- 
stand the alteration of the AEoi term we should take 
into consideration other components of the orbital inter- 
actions, such as the orbital polarization and the interlayer 
interaction of the tt-tt* orbitals. 

Orbital polarization reflects a mixing of the occu- 
pied/virtual orbitals in one fragment due to the pres- 
ence of another fragment, i.e., each valence electron of 
one fragment entering the electron space of electrons of 
other fragment polarizes its orbitals. The polarization 
effect is caused by the repulsion arising between inter- 
acting electrons [26]. Analyzing the orbital formation 
after perturbation of the fragment's orbitals, a discrep- 
ancy in the product orbitals for the AB, AA and AA' 
stacking patterns was detected as demonstrated by the 
scheme in Fig. [2] Let us consider the formation of the 
orbitals generating the HOMO-LUMO gap in the joint 
system which would give the most contribution into the 
interlayer interaction of the tt-tt* orbitals. In formation 
of HOMO orbital (LUMO) of the joint system two degen- 
erate fragment's orbitals HOFO participates (two LUFO 
orbitals for the LUMO formation). Perturbation of those 
degenerate HOFO orbitals (overall four in two fragments) 
creates the four molecular orbitals in the bilayer system 
(HOMO, HOMO-1, HOMO-2 and HOMO-3). Simi- 
larly, the LUFO orbitals perturb in the valence band so 
that two LUFOs are taken from each fragment and their 
perturbation leads to formation of four LUMO orbitals 
in the stacked system (LUMO, LUMO+1, LUM0-K2, 
LUMO-l-3). For the AA stacking, two pairs of product 
orbitals possess an identical orbital energy, i.e., HOMO 
and HOMO-1 (LUMO and LUMO+1) are degenerate, 
the same for HOMO-2 and HOMO-3 (LUM0-f2 and 
LUMO-l-3). However, already for the AB case the con- 
duction band is limited by a single HOMO orbital being 
a product of the perturbation of all four fragment's or- 
bitals while the rest of the generated orbitals are shifted 
deeper into the conduction band where two of them still 



would remain degenerate. Mixing of the LUFOs for the 
AB stacking stays similar to those for the AA case, i.e., 
two pairs of the degenerate orbitals are formed. For the 
AA' stacking all four product orbitals HOMOs (LUMOs) 
are separated by the energy gap. 



The spatial orbital overlap J.^ 






regardless of the 



observed peculiarities of orbital mixing as lattice ar- 
rangement is changed, is being affected mostly by re- 
arrangements of the TT clouds from their superposition 
in AA stacking. However, an analysis of interaction 
of the occupied/unoccupied orbitals between fragments 
J^f^ has shown a distinct behavior. We observed barely 
present overlap between those orbitals in the AA pattern 
iJi^j^^= 10~^ eV), while it appears for the AB stacking 
to be 0.169 eV and increases even further up to 0.191 eV 
for the AA' stacking. Such progress is consistent with 
improvement of the orbital interactions AEoi and with 
enhancement of the attractive interactions of the electro- 
static nature (see AVei in Table |I| as stacking pattern 
is changed from AA to AB. Both these terms, AEoi and 
AVei, lowers the bonding energy AE'^ and their contri- 
bution compensate the growing Pauli repulsion between 
fragments. 

Therefore, the efficiency of the acceptor-donor interac- 
tions, defined by AEoi, is found to be several times (at 
least five) weaker for the AA stacking in comparison to 
that for the AB. This behavior offers an interpretation 
of the interlayer decoupling observed in the experiments 
for the AA stacking in twisted bilayer graphene [iHal, 
against the efficient coupling known for the AB stacking. 
According to our findings the interlayer decoupling in 
AA staking must be caused by suppression of the orbital 
overlap of the occupied, tt, and unoccupied, tt*, orbitals 
between the fragments, i.e. inter-band (HOMO-LUMO) 
interaction, while variation of the tt-tt overlap is found 
to have a insignificant affect on the interaction between 
fragments. This conclusion is, in fact, in contrast to the 
commonly accepted opinion that the electronic coupling 
between the stacked graphene layers, which is capable 
to change linear Dirac cones dispersion to parabolic one, 
originates from intra-band interlayer interactions [7, 25] . 
We should also note that along with the intra- and inter- 
band interactions the effect of orbital polarization is one 
of the significant contribution into the bonding interac- 
tions. 



III. INTERLAYER DISTANCE 

The observation that the interlayer repulsion charac- 
terized by AEp is being stronger for the AB stacking 
in comparison to that for AA was beyond our expecta- 
tions because the Pauli repulsion is presumed to domi- 
nate when one structure of closed electron shell is placed 
exactly on top of the other (such as the coordinates in 
x — y directions coincide as presented in Fig. [U (a)). The 
Pauli repulsion has the exponential dependence on the 
separation distance and an increase in deq for the AA 
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FIG. 2: The energy diagram demonstrating the formation of the molecular orbitals (HOMOs and LUMOs) as the fragment's 
orbitals overlap (HOFO and LUFO) for different arrangements of graphene lattices in the case when the fragments are separated 
by the equilibrium distance dec- 



ease must cause a significant suppression of the repul- 
sive forces. To understand this, we considered the de- 
viation of the bonding energy and its components with 
the change in the interlayer distance for the AA stack- 
ing [Fig. [T] (b)]. The negative sign is for the attractive 
interactions. The electrostatic interactions AVei and the 
orbital interactions AEoi are both attractive in nature 
and their magnitudes are reduced with increasing dis- 
tance between the fragments (AVei can have a positive 
magnitude only for short interlayer distance when the 
nuclear repulsion dominates over other attractive terms) . 
The contribution of the repulsive forces collected within 
AEp into the total bonding energy displays its domina- 
tion over the attractive terms only for a short interlayer 
distance, so the AE"" remains positive (repulsive) up to 
a distance of 3.2 A. The bonding energy AE^ reaches its 
minimum at deq '-3.67 A {AE° ~ -1.47 eV) which is still 
energetically far from the ground state of AB stacking, 
where AE° = -1.93 eV. 

Among all the interaction terms the Pauli repulsion 
and the bonding interactions are points of particular in- 
terest. According to the Pauli principle the valence elec- 
tron from one fragment is not supposed to penetrate the 
closed valence shell of the other fragment because the re- 
pulsive forces expel the charges from the overlap region. 
The main component within AEp which contribute to the 
repulsion effect comes from the kinetic energy while the 
potential energy part is attractive [21|. It is noticeable 
in Fig. [1] (b) that as the fragments being separated at a 
distance beyond the value of 4.0 A the Pauli repulsion 
becomes negligible. However, for the interlayer distance 
3.4 A which is typical for the AB stacking, the Pauli re- 
pulsion is still strong as A_Ep=4.33 eV while for the AB 
stacking pattern it was 3.38 eV. Therefore, if we compare 
the AA and AB stacking for the same interlayer distance 
3.4 A, the Pauli repulsion is larger by almost 1.0 eV for 
the AA stacking, in agreement with our expectations, as 



stated above. 

The orbital interaction energy AEoi belongs to the at- 
tractive forces and its value reflects the efficiency of the 
donor-acceptor charge transfer between fragments which 
is controlled by the interlayer orbital overlap together 
with the Pauli repulsion and orbital polarization. Sup- 
pression of the Pauli repulsion with growing interlayer 
distance tends to increase the orbital interaction energy 
AEoi- However, a simultaneous reduction of the spatial 
orbital overlap generally leads to diminishing of AEoi, 
i.e., to a reduction of the charge transfer between the 
fragments. As a result, AEoi being strongly attractive 
(negative sign) at short distances almost vanishes as the 
distance reaches the value of d ~ 3.3 A while after d ~3.9 
A it turns repulsive with a positive sign. It was noticed 
that the composition of the molecular orbitals near the 
HOMO-LUMO gap is not changed with distance as it is 
shown in Fig. [S] For example, if the HOMO was formed 
by the interaction of selective HOFO orbitals provided 
by each fragment it remains of the same composition on 
a whole range of the distance while just become shifted 
in energy. Therefore, we can conclude that the orbital 
polarization term brings no contribution in deviation of 
AEoi with distance. 

To separate the orbital interactions from other forces, 
we calculated the spatial overlap integral between the 
fragments J^-"-^ and j"~^- The degradation of j"f" 
upon increase of the interlayer distance along with the 
size of the HOMO-LUMO gap [Egap) are presented in 
Tableini Because a rise in interlayer separation d induces 
a suppression of the orbital overlap, in particular of the 
overlap matrix 5'ij-, the charge transfer integral J^f^ 
is also being reduced. The same gradual reduction is 
observed for the Jj ~ overlap (from 5xl0~^ to 6xl0~® 
eV). 

As one increases the interlayer distance d the HOMO- 
LUMO gap grows in contrast to the diminishing Jij and 



TABLE II: The spatial overlap integral j"~" and HOMO- 
LUMO gap AEgap calculated for the AA-stacking pattern as 
the flake separation d gradually increases. All values are in 
eV. 



d 


2.94 


3.f4 


3.34 


3.40 


3.54 


3.74 


3.94 


4.f4 


jH-H 

AEgap 


-f.2f 
0.45 


-0.91 
0.84 


-0.68 
f.f3 


-0.63 
f.20 


-0.5f 
f.34 


-0.38 
f.50 


-0.28 
f.6I 


-0.2f 
f.69 



its enhancement is directly connected to the orbital in- 
teraction between the fragments. We presented in Fig. [3] 
the energy diagram for the energetics of the n and tt* 
orbitals near the HOMO-LUMO gap for the case of sep- 
arated fragments and their orbital splitting/mixing af- 
ter perturbation. To find the orbital energy change by 
the fragment interaction we used the expression derived 
within the Hiickel approximation for the description of 
splitting of the tt orbitals belonging to different fragments 
after inclusion of the interactions 



Si « eo + Hi2 - (eo + Hi2)Si 



E2 ~ So — H12 + (eo — H 12)812 



(1) 



(2) 



where eo and Ei^2 are the molecular 7r-orbital energies 
before (for identical graphene flakes ei = 62 = eg) and 
after perturbation, respectively. S12 is the orbital overlap 
between the fragments and H12 is the intrinsic interac- 
tion integral that is a combination of the energy terms 
responsible for electron-electron interactions and partic- 
ularly the contribution of its repulsive part to the orbital 
energies £^1(2). 

As fragments are brought together to a distance of 2.94 
A (see Fig. [3]) which is much shorter than the equilibrium 
separation {df,q — 3.67 A), the repulsive force dominates 
over the attractive interaction. Therefore, since the elec- 
trons are expelled by repulsion from the overlap region, 
within the Hiickel approximation the overlap matrix 5*12 
is treated as zero or would possess a negative value so 
that for simplicity, we can ignore the contribution from 
(eo ± Hi2)Si2 to the orbital energy i?i(2)- Besides the 
effect of repulsion on (eo ± Hi2)Si2 term, the strong re- 
pulsion between the fragments causes the increment of 
the H12 and therefore, closer the fragments are to each 
other the larger the splitting of their tt orbitals. The sig- 
nificant TT — TT splitting leads to shift of the LUMO and 
HOMO orbitals close to each other thus causing a sup- 
pression of the HOMO-LUMO gap. However, for a large 
interlayer distance of 4.14 A, the value of H12 dimin- 
ishes because of the suppression of the Pauli repulsion, 
while the overlap matrix would have positive values. The 
charge exchange between the fragments is allowed which 
decreases the splitting \Ei— E2\ for the tt orbitals and in 
turn enlarges the HOMO-LUMO gap {/AEgap). 

In this section we considered the basics of the orbital 
interactions between two graphene fragments being un- 
der the control of the separation distance. It was gath- 
ered that there are two main components affecting the 



efficiency of the orbital interaction energy by varying the 
distance: the Pauli repulsion and the spatial orbital over- 
lap, while no contribution from the orbital polarization 
was observed. The spatial orbital overlap decreases with 
increasing distance and so is the Pauli repulsion. These 
changes have the opposite influence on the orbital inter- 
action energy AEoi , but since the decrease of the spatial 
orbital overlap is faster than the Pauli repulsion, AEoi is 
generally reduced with induced flake separation. 



IV. FLAKE ROTATION 

For two graphene layers stacked together, the rota- 
tional misorientation creates the Moire pattern which is 
a periodic pattern manifesting itself through the spots 
where the superposition of the lattices is preserved, i.e., 
the AA stacking |11-13]. The rest of the surface (the 
inter-spot regions) which are of a larger amount of the 
surface, possesses the stacking order similar or close to 
that for the AB arrangement (see Ref. [23| for images of 
various rotation structures). As the rotation angle in- 
creases, the percentage of AA-spotted areas grows while 
the size of the spots and the inter-spot areas shrink. Al- 
ready for the angle above 20° the spots with ideal AA 
stacking vanish completely, leaving the mixed and weakly 
defined interlayer lattice order. Therefore, the experi- 
mentally observed lattice misorientation induced by ro- 
tation angle of 3-5° (see Refs.[l]-Q) is characterized by 
the well defined the AA-spots and of a high coverage of 
the AB stacking. 

Graphene flakes of finite size are considered in this 
paper as a model system to represent the spots char- 
acterized by the lattice superposition (the AA-spots) in 
misoriented bilayer graphene. Therefore, by rotating one 
flake with respect to the other (the rotation axis is placed 
at the flake center) we basically recreate the modification 
of the shape of the AA-spots affected by the rotation an- 
gle. The main disadvantage of describing misoriented 
bilayer graphene by finite size flakes is to underestimate 
the contribution of the AB stacking areas to the elec- 
tronic properties of the bilayer system. In the model sys- 
tem of rotated flakes the contribution of the inter-sport 
regions with AB stacking would depend on the flake size 
and therefore for small flakes would be negligible. 

In a system of the AA stacked graphene fragments, 
the interacting tt orbitals between the fragments are per- 
fectly orthogonal. The layer rotation observed in the ex- 
periments [J-Q breaks that orthogonality, thereby mod- 
ifying the balance of the attractive and repulsive forces 
between the layers. We simulated the effect of misorien- 
tation in the model system of two graphene flakes pre- 
sented in Fig. [T] (a) with the rotation axis placed at the 
flake center. In a system of two flakes stacked in the AA 
pattern and separated by a distance of d=3.4 A which is 
the equilibrium distance for the AB-stacking, we rotate 
one flake relative to the other for further elaboration of 
the interaction parameters. Our simulation results for 
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FIG. 3: Energetics of the n and n* orbitals defining the HOMO-LUMO gap in the separated fragments and their sphtting as 
fragments being stacked in the AA lattice arrangement. 



TABLE III: The interlayer interactions in the system of two 
twisted flakes separated by a distance d—3A A. For the spatial 
overlap integral J^~^ , its modulus has been considered to 
avoid confusions of the sign change. /S.Egap is the HOMO- 
LUMO gap in the system of two stacked flakes. For a better 
comparison we repeated the results for the graphene flakes 
stacked in the AB pattern at the equilibrium distance deq—3.4 
A. All values are in eV 
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the bonding energy /S.E^ between layers and its compo- 
nents, overlap integrals such as Jj " and Jj a' , and the 
HOMO-LUMO gap (AE'gap) are displayed in Table [ml 
The resuhs indicate that the size of the HOMO-LUMO 
gap {/S.Egap) grows as the bonding interaction between 
layers is improved with rotation. 

When two lattices of different flakes are superposed in 
the space, the electronic clouds of their tt orbitals are 
also superposed giving the maximum magnitude of the 
overlap integral Jj " which is being suppressed with the 



flake rotation because of misorientation of those tt clouds. 
Therefore, the spatial overlap integral Jj " reaches its 
minimum as the rotation angle reaches the value of ^ ~ 
18° which is being a result of significant disarrangement 
of the TT orbitals within the overlapping region from their 
superposed position. Another overlap integral account- 
ing for inter-band interactions, J/^-~^ (between HOMO 
and LUMO orbitals belonging to the different fragments), 
has shown the opposite behavior to that of J^~^ , i.e., 
with its minimum for the superposed case while grow- 
ing with the flake rotation. However, even when J^~^ 
reaches the maximum bX 8 — 24°-30°, its magnitude is 
still much lower than that found for the AB-stacking. 
For the angle in the range oi = 24°-30°, both integrals 
{Ji^r^ and Jj^f^) deviate insignificantly. With break- 
ing of the orthogonality of the tt orbitals as the flake is 
rotated, the Pauli repulsion is also being suppressed by ~ 
0.5 eV when its reaches its minimum at ~ 18°. The re- 
duced value of 3.831 eV for 9 ~ 18° becomes much closer 
to that for the AB stacking (A^p=3.388 eV). 

Moreover, the orbital interaction energy AEoi grows 
with flake rotation as the inter-band interactions reflected 
by the J.^j^^ improve {J^r^ increases up to several or- 



ders of magnitude) along with the fast reduction of the 
interlayer repulsion. The coincidence of the minimum 
of A_Eoi with the minimum value of the Pauli repulsion 
achieved for the 9 ~ 18° is a clear evidence that the ef- 
ficiency of the orbital interactions is being under direct 
control of the repulsive forces. 

Therefore, because the layer rotation significantly sup- 
presses the interlayer repulsion which in turn improves 
the orbital interaction, regardless of the deviation of 
the overlap integrals, the total bonding energy AiJ° 
is lowered with rotation and its magnitude reaches its 
minimum also at ~ 18°. In fact, the magnitude of 
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the total bonding energy found for the angle 6 ~ 18° 
(A£;°:=-1.837 eV) is comparable to that for the AB 
stacking (A£;°=-1.931 eV in Table U). Above the rota- 
tion angle 30°, some fluctuations for all the terms occur 
while closer to angle 60° for which the conditions for the 
lattice superposition between two layers reappear, i.e., all 
the terms have the same values as for the angle 0°. Ba- 
sically, the dependence in the range of the rotation angle 
from 30° to 60° is displayed in reverse order to that from 
0° to 30°. 

Our main conclusion is that the repulsion appeared as 
result of the interaction of two systems of closed electron 
shells whose lattices are superposed, is the central force 
controlling the efficiency of the interlayer orbital inter- 
actions. As the flake rotation breaks the lattice super- 
position, the suppression of the repulsion induces an im- 
provement of orbital interactions (such as the electronic 
coupling) thereby lowering the bonding energy lS.E^ . Be- 
cause the orbital interaction between the flakes depends 
on the flake rotation, the equilibrium distance deq also 
fluctuates with rotation. The modification of the equi- 
librium distance d^q follows the dependence observed for 
the bonding energy AE^ which is controlled by the Pauli 
repulsion. The maximum equilibrium distance (ieij=3.67 
A is obtained for the AA stacking which is suppressed 
down to deg=3.43 A as the flakes are at an angle 18° 
that brings the system to the lowest energy state achiev- 
able with rotation. If the rotational angle grows further 
beyond 18°, the magnitude of deq enhances again and for 
e ~ 30° its value is 3.54 A. 



V. DISCUSSION AND SUMMARY 

Two graphene layers stacked in the AA pattern which 
is characterized by the lattice superposition between the 
layers is the most unstable configuration in the bilayer 
geometry. The instability appears as a result of strong in- 
terlayer repulsion induced by the interaction of the filled 
orthogonal tt orbitals within the overlap region. There- 
fore, a disruption of the lattice superposition lowers the 
total bonding energy and therefore, leads to an enhance- 
ment of the system stability. The AB stacking is the 
most successful scheme to suppress the Pauli repulsion 
because it induces the maximum mis-orientation in the 
interlayer lattice order and therefore, the AB stacking 
appear to be the ground state of the system character- 
ized by the strongest bonding interactions between the 
layers. However, an alternative way to induce the lattice 
mis-orientation from that superposed in the AA stacking 
and thereby, to transfer the system to the lower energy 
state is the layer rotation. It should be noted that any 
modification of the lattice order different from the AB- 
stacking would be metastable (such as AA') because AB 
stacking is the ground state of the system. 

For adjacent graphene flakes of finite size, the rotation 
of one flake relative to the other induces a fast reduc- 
tion of the repulsive part (Ai?p) and an increase of the 



attractive forces (orbital interaction energy lS.Eoi) such 
that both these tendencies lower the total bonding en- 
ergy between the layers /^E^ . The lowest bonding en- 
ergy Ai?° is achieved for the rotation angle of 6* ~ 18°. 
This state is still metastable but with the lowest value of 
the bonding energy among all the rotation angles, and 
its value correlates with that for the AB-stacking be- 
ing the ground state of bilayer graphene. However, as 
it was already noted above, the description of the adja- 
cent graphene layers by the model system of finite flakes 
has crucial disadvantage caused by the underestimation 
of the inter-spot areas of the AB-stacking into the in- 
terlayer repulsion. As we switch to the twisted bilayer 
graphene of the infinite size, the contribution of large ar- 
eas of AB stacking, which was largely neglected in the 
fiake system, should be considered. For a small rotation 
angle the percentage of the AA spotted areas is large 
which decreases with angle enhancement. Thus, for an 
angle altered from 10° to 12°, a decrease of 2.5 % of AA 
staking is observed [23| . Obviously, since the larger inter- 
spot areas of AB-stacking is observed for small rotation 
angles 2°-5° (see Ref.[23| for the images of various rota- 
tion structures), we would expect that the rotation angle 
of much smaller magnitude than that for the finite sys- 
tems might be required {9 < 18°) to bring the stacked 
graphene layers to the metastable state with the lowest 
energy. 

However, additionally to rotation there is another way 
to make the system of the fiakes stacked in the AA ar- 
rangement more stable which is to raise the interlayer 
distance when the bonding energy is lowered again due 
to suppression of the Pauli repulsion between the layers 
(see Fig. [1] (b)). In case of misoriented bilayer graphene 
exhibiting the Moire pattern the distribution of the re- 
pulsive forces would be non-uniform as the lattice order 
is not the same in different areas, i.e. a maximum force 
pushing apart two lattices would originate at the AA- 
spots of the lattice superposition. Moreover, another 
interesting distinction between the system of adjacent 
fiakes and bilayer graphene of infinite size is alteration of 
its rigidity. Recalling that the free standing graphene is 
subjected to rippling of its surface [28| , the lower rigidity 
of the graphene layers than that of flakes is anticipated. 
Therefore, we expect that a strong Pauli repulsion which 
is pronounced locally at the center of the AA-spots might 
not able to modify the interlayer separation throughout 
the whole system because of large areas of AB-stacking, 
but would rather induce a local lattice distortion forming 
the bump on the surface with its highest point at the cen- 
ter of the AA-spot [IS]. To simulate this effect, the fiake 
of larger size containing the bigger areas of AB stack- 
ing have been examined and already for that system the 
generation of the bump as high as 0.2 A was observed. 
However, the bump's height may be enhanced for an in- 
finite system due to better efficiency of the attractive 
interactions in the inter-spot areas and lower rigidity of 
the layers. In fact, the appearance of bumps can explain 
the brightening of the AA-spots observed in the STM 



images of the twisted bilayer graphene [ll-Q • 

The final point we wish to make is about the inter- 
layer coupling which is found to be a function of the ro- 
tation angle and interlayer distance. As fragments being 
stacked, there are two types of tt interactions occurs, such 
as perturbation of the occupied/occupied orbitals (tt-tt 
interaction) and interaction of the occupied/unoccupied 
orbitals (tt-tt* interaction). According to the theoreti- 
cal models developed to describe the behavior of the tt 
bands in bilayer graphene, the modification of the linear 
dispersion of the Dirac cones to a parabolic one has been 
simulated by inclusion only of the tt-tt orbital interaction 



[3,[23- However, we found that although the intra-bands 
interactions plays an important role, but particularly the 
inter-bands part must be introduced in the model to ac- 
count for the decoupling effect arising in the AA stacked 
graphene layers. The orbital interaction energy AEoi is 
therefore suppressed at least by five times as the lattice 
arrangement was changed from AB to AA due to van- 
ishing of the TT-TT* interactions between fragments that 
is reflected by a drastic reduction of the spatial overlap 
integral j"~^- 
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